//==================================================================================================
/*
  EVE - Expressive Vector Engine
  Copyright : EVE Project Contributors
  SPDX-License-Identifier: BSL-1.0
*/
//==================================================================================================
#include "test.hpp"

#include <eve/module/bessel.hpp>

TTS_CASE_TPL ( "Check eve::cyl_bessel_k over real, positive floating orders"
             , eve::test::scalar::ieee_reals
              )
<typename T>(tts::type<T>)
{
  auto constexpr N = 11;
  using a_t  = std::array<T, 8 >;
  using aN_t = std::array<a_t, N >;
  a_t re{ 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00,     3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
  aN_t resN8;
  //res are taken from octave 9.2.0 besselk outputs
  int p=0;
  resN8[p] = a_t{1.5091129245821360e+00,4.3843063344153255e-01,1.6992206277161609e-01,7.1612291089370450e-02,3.1444799766003367e-02,1.4149839171199918e-02,6.4720581217078289e-03,2.9950130446850290e-03,};
  ++p; //1  1.333e+00
  resN8[p] = a_t{5.3402033678005392e+00,7.8676215106523129e-01,2.5048339894814947e-01,9.5843815923198697e-02,3.9721125205585531e-02,1.7197900351026186e-02,7.6521425170157599e-03,3.4684514504195935e-03,};
  ++p; //2  2.333e+00
  resN8[p] = a_t{4.8977587305031470e+01,2.5364630362821523e+00,5.6283719837655766e-01,1.7810541989292461e-01,6.5613509620270502e-02,2.6218541171920049e-02,1.1006661094754209e-02,4.7737060961822566e-03,};
  ++p; //3  3.333e+00
  resN8[p] = a_t{7.6721378366829049e+02,1.2623589653715277e+01,1.7955266886092860e+00,4.4215991015944112e-01,1.3849415044040136e-01,4.9396108807770102e-02,1.9066457726390498e-02,7.7525466649421317e-03,};
  ++p; //4  4.333e+00
  resN8[p] = a_t{1.7098172779933688e+04,8.6693727394383842e+01,7.6041183301776618e+00,1.4063273925580384e+00,3.6345039228780018e-01,1.1287838118555174e-01,3.9253265133851245e-02,1.4712868487133704e-02,};
  ++p; //5  5.333e+00
  resN8[p] = a_t{4.9471442742619739e+05,7.6396922707170847e+02,4.0561620136573829e+01,5.5205643832856914e+00,1.1545920213525309e+00,3.0683803080990552e-01,9.4665338724918849e-02,3.2273994143498304e-02,};
  ++p; //6  6.333e+00
  resN8[p] = a_t{1.7606944481266961e+07,8.2356988161592744e+03,2.6210840154005263e+02,2.5942169096050002e+01,4.3362401431782294e+00,9.7417811679230371e-01,2.6364517914847374e-01,8.0915933396873807e-02,};
  ++p; //7  7.333e+00
  resN8[p] = a_t{7.4389903696980953e+08,1.0508282089842251e+05,1.9935261806310832e+03,1.4243756794577183e+02,1.8872562498854979e+01,3.5540984201175836e+00,8.3677769484654896e-01,2.2937690882819090e-01,};
  ++p; //8  8.333e+00
  resN8[p] = a_t{3.6386004307449890e+10,1.5494504053263613e+06,1.7461157803063150e+04,8.9639397320910143e+02,9.3625783148513790e+01,1.4691750966368954e+01,2.9909206290187123e+00,7.2787644547638686e-01,};
  ++p; //9  9.333e+00
  resN8[p] = a_t{2.0221885827841873e+12,2.5929256243004445e+07,1.7318134777928935e+05,6.3673957152312005e+03,5.2223698802849924e+02,6.7991602658577889e+01,1.1914261506026971e+01,2.5623142340730207e+00,};
  ++p; //10  1.033e+01
  resN8[p] = a_t{1.2586145337754586e+14,4.8556223360807604e+08,1.9190602706736906e+06,5.0420582869451791e+04,3.2382786142878654e+03,3.4868558858745320e+02,5.2413042431797273e+01,9.9259275421487683e+00,};
  using eve::spherical;
  auto v0 = T(1.0/3.0);
  auto h =  eve::half(eve::as<T>());
  for(int j=0; j <= 7; ++j)
  {
    auto c = re[j];
    auto fac = eve::sqrt(eve::pio_2(eve::as(c))*eve::rec(c));
    for(int i=0; i < N; ++i)
    {
      auto v =  v0+i;
      auto res = resN8[i][j];
      TTS_RELATIVE_EQUAL(eve::bessel_k(v, c), res, tts::prec<T>());
      TTS_RELATIVE_EQUAL(eve::bessel_k[spherical](v, c), eve::bessel_k(v+h, c)*fac, tts::prec<T>());
    }
  }
};

TTS_CASE_TPL ( "Check eve::cyl_bessel_k over real, negative floating orders"
             , eve::test::scalar::ieee_reals
              )
<typename T>(tts::type<T>)
{
  auto constexpr N = 11;
  using a_t  = std::array<T, 8 >;
  using aN_t = std::array<a_t, N >;
  a_t re{ 2.9999999999999982e-01, 1.0000000000000000e+00, 1.7000000000000002e+00, 2.3999999999999995e+00,     3.0999999999999988e+00, 3.7999999999999998e+00, 4.4999999999999991e+00, 5.2000000000000002e+00};
  aN_t resN8;
  //res are taken from octave 9.2.0 besselk outputs
  int p=0;
  resN8[p] = a_t{1.5091129245821360e+00,4.3843063344153255e-01,1.6992206277161609e-01,7.1612291089370450e-02,3.1444799766003367e-02,1.4149839171199918e-02,6.4720581217078289e-03,2.9950130446850290e-03,};
  ++p; //-1 -1.333e+00
  resN8[p] = a_t{5.3402033678005392e+00,7.8676215106523129e-01,2.5048339894814947e-01,9.5843815923198697e-02,3.9721125205585531e-02,1.7197900351026186e-02,7.6521425170157599e-03,3.4684514504195935e-03,};
  ++p; //-2 -2.333e+00
  resN8[p] = a_t{4.8977587305031470e+01,2.5364630362821523e+00,5.6283719837655766e-01,1.7810541989292461e-01,6.5613509620270502e-02,2.6218541171920049e-02,1.1006661094754209e-02,4.7737060961822566e-03,};
  ++p; //-3 -3.333e+00
  resN8[p] = a_t{7.6721378366829049e+02,1.2623589653715277e+01,1.7955266886092860e+00,4.4215991015944112e-01,1.3849415044040136e-01,4.9396108807770102e-02,1.9066457726390498e-02,7.7525466649421317e-03,};
  ++p; //-4 -4.333e+00
  resN8[p] = a_t{1.7098172779933688e+04,8.6693727394383842e+01,7.6041183301776618e+00,1.4063273925580384e+00,3.6345039228780018e-01,1.1287838118555174e-01,3.9253265133851245e-02,1.4712868487133704e-02,};
  ++p; //-5 -5.333e+00
  resN8[p] = a_t{4.9471442742619739e+05,7.6396922707170847e+02,4.0561620136573829e+01,5.5205643832856914e+00,1.1545920213525309e+00,3.0683803080990552e-01,9.4665338724918849e-02,3.2273994143498304e-02,};
  ++p; //-6 -6.333e+00
  resN8[p] = a_t{1.7606944481266961e+07,8.2356988161592744e+03,2.6210840154005263e+02,2.5942169096050002e+01,4.3362401431782294e+00,9.7417811679230371e-01,2.6364517914847374e-01,8.0915933396873807e-02,};
  ++p; //-7 -7.333e+00
  resN8[p] = a_t{7.4389903696980953e+08,1.0508282089842251e+05,1.9935261806310832e+03,1.4243756794577183e+02,1.8872562498854979e+01,3.5540984201175836e+00,8.3677769484654896e-01,2.2937690882819090e-01,};
  ++p; //-8 -8.333e+00
  resN8[p] = a_t{3.6386004307449890e+10,1.5494504053263613e+06,1.7461157803063150e+04,8.9639397320910143e+02,9.3625783148513790e+01,1.4691750966368954e+01,2.9909206290187123e+00,7.2787644547638686e-01,};
  ++p; //-9 -9.333e+00
  resN8[p] = a_t{2.0221885827841873e+12,2.5929256243004445e+07,1.7318134777928935e+05,6.3673957152312005e+03,5.2223698802849924e+02,6.7991602658577889e+01,1.1914261506026971e+01,2.5623142340730207e+00,};
  ++p; //-10 -1.033e+01
  resN8[p] = a_t{1.2586145337754586e+14,4.8556223360807604e+08,1.9190602706736906e+06,5.0420582869451791e+04,3.2382786142878654e+03,3.4868558858745320e+02,5.2413042431797273e+01,9.9259275421487683e+00,};
  using eve::spherical;
  auto v0 = T(-1.0/3.0);
  auto h =  eve::half(eve::as<T>());
  for(int j=0; j <= 7; ++j)
  {
    auto c = re[j];
    auto fac = eve::sqrt(eve::pio_2(eve::as(c))*eve::rec(c));
    for(int i=0; i < N; ++i)
    {
      auto v =  v0-i;
      auto res = resN8[i][j];
      TTS_RELATIVE_EQUAL(eve::bessel_k(v, c), res, tts::prec<T>());
      TTS_RELATIVE_EQUAL(eve::bessel_k[spherical](v, c), eve::bessel_k(v+h, c)*fac, tts::prec<T>());
    }
  }
};
